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Superconductivity in the two-dimensional Hubbard model? 

D. Baeriswyl*, D. Eichenberger, and B. Gut 

Department of Physics, University of Fribourg, Chemin du Musee 3, CH-1700 Fribourg, Switzerland 

Received zzz, revised zzz, accepted zzz 
Published online zzz 

PACS 71.10.Fd, 74.20.Mn 

A refined variational wave function for the two-dimensional repulsive Hubbard model is studied numeri- 
cally, with the aim of approaching the difficult crossover regime of intermediate values of U. The issue of a 
superconducting ground state with d-wave symmetry is investigated for an average electron density n ~ 0.8 
and for U = 8t. Due to finite-size effects a clear-cut answer to this fundamental question has not yet been 
reached. 
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1 Introduction Strong correlations between electrons play a fundamental role in the cuprates. In fact, 
there is ample evidence for a crossover from a doped Mott insulator to a correlated Fermi liquid, as holes 
are introduced into the G1O2 planes. Two limiting regimes can be described in simple terms, an antifer- 
romagnetic Mott insulator at half filling (x = 0) and a rather conventional Fermi liquid at high doping 
levels (x > 0.25). Many difficult problems remain to be solved, especially concerning the nature of the 
underdoped region (the so-called pseudogap phase), the normal phase around optimal doping ("marginal 
Fermi liquid") and the "physical mechanism" leading to the superconducting phase for 0.05 < x < 0.25. 
In this note we discuss two specific issues: on the one hand, the nature of the crossover between the under- 
doped and the overdoped regimes, on the other hand, the possibility of superconductivity originating from 
purely repulsive interactions. We use variational wave functions for the one-band Hubbard model in two 
dimensions, keeping in mind that both the model and the method do provide insight, but at the same time 
are not sufficient for developing a complete theory of the layered cuprates. 

The Hubbard Hamiltonian H = —tT + UD is composed of two terms with conflicting tendencies, the 
hopping term 

f = E (i) 

where cL creates an electron at site i with spin er, and the number of doubly occupied sites 

D = ^2,n^mi, (2) 

i 

where m a = c\ a Ci a . The notation (i, j) means that hopping is restricted to neighboring sites. 

If the on-site repulsion U vanishes, the electrons occupy delocalized Bloch orbitals and set up a filled 
Fermi sea in the ground state \ipo), given by 

n>o>= n c Lio). o) 

k,cr,e k <e J! . 
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In contrast, for a vanishing hopping amplitude t (and U > 0, or for finite t and U — > oo ), the ground 
state | ^00} for N = 2M electrons is an arbitrary superposition of configurations where the electrons are 
localized on lattice sites without any double occupancy, 

IV'oo) = Yl ^ u " ' ' ^h, ■ ■ ■ Jm) ct i)T • ■ • 4 M ,Ai,l ■ ■ ■ c Lj°) ' ik ^ h ■ (4) 
ii, ■ ■ ■ , hi 
h , • • ■ . hi 

For finite values of U and t the delocalizing tendency of t and the localizing tendency of U compete. One 
expects a crossover to occur around some critical value (U/t) c between a ground state linked to \tpo) for 
U < U c and one linked to \tpoo) for U > U c - This picture is captured by two types of variational wave 
functions, the Gutzwiller ansatz H] 

\ V ) G ^e^ B \ Vo ) (5) 
and its counterpart f2| 

\v)b = e~ hf l^oo) , (6) 
where g and h are variational parameters. Depending on the relative values of the minima of 

(Wg) r „ 

Waive) {Wb\Wb) 
one or the other of the two ground states is favored. 



2 Mott transition and localization For the case of a half-filled band, n = N/ L = 1 (L is the total 
number of sites), the two variational states defined above can be distinguished by their sensitivity with 
respect to changes in boundary conditions. Using an early argument of Kohn J5), it can be readily shown 
that \ipc) is metallic (the Drude weight is finite) and \4>b) is insulating (vanishing Drude weight) [4|. Thus 
the delocalization-localization crossover manifests itself in this approach as a sharp (Mott) metal-insulator 
transition at a critical value U c , which turns out to be of the order of the bandwidth. The distinction 
between the two regimes is blurred if we allow for antiferromagnetic ordering, which is expected to set in 
at arbitrarily small U for the square lattice. In this case, a possible signature would be the size of the local 
moment, which is small in the "delocalized" spin-density-wave regime at small U, increases strongly in 
the crossover regime and saturates for U — > 00. 

The insulating cuprates can be well understood in terms of the antiferromagnetic Heisenberg model, the 
strong coupling limit of the half-filled Hubbard model, provided that ring exchange around the plaquettes 
is also included. The presence of this term indicates that U is not much larger than the bandwidth. In fact, 
the analysis of the spin waves observed in neutron scattering experiments yields a value of U of the order 
of 12t (twice the bandwidth) [5 1. According to our variational analysis this puts these materials inside the 
localized regime, but not so far that double occupancy can be simply discarded. At half filling, cuprates 
can thus be considered as "weakly localized" Mott insulators. The good agreement between theoretical 
predictions made on the basis of the (two-dimensional) Hubbard model and the observed magnetic exitation 
spectrum in the layered cuprates does not imply that long-range Coulomb interactions are absent in these 
materials; these are just not seen at half filling. The addition of non-local Coulomb terms would essentially 
lead to slightly different exchange constants. 

An interesting and largely unexplored issue is the location of the crossover regime (as a function of U) 
away from half filling. The variational procedure should also be useful in this limit, although the presence 
of holes poses an additional challenge. The distinction between the localized and delocalized regimes is 
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less clear-cut for n < 1 than it is at half filling because the Drude weight is expected to remain finite 
due to the motion of holes. A further complication arises from the fact that the Hubbard Hamiltonian 
represents an effective model, in which many degrees of freedom are subsumed in some average way by 
the parameters. Quite generally, we expect that screening will become more and more efficient as doping 
is increased and thus will lead to a reduced value of U. The appearance of a conventional Fermi liquid 
above a hole concentration of about 25% agrees with such a picture. 



3 rf-wave superconductivity The celebrated BCS theory of superconductivity is based on an effective 
attraction between electrons, caused by phonon exchange. Therefore it is not surprising that the application 
of the BCS mean-field approximation to the repulsive Hubbard model does not yield any energy gain due 
to a pairing instability. To see this, we decompose the on-site term, 

(n iT n a ) « (nif)(nn) + (c] T c iL )(c iT ct i ) + {c^cj^ia^) , (8) 

and notice that in the absence of antiferromagnetism the first term is just \n 2 , while the second term 
vanishes. The third term, responsible for a superconducting instability in BCS theory, is nonnegative, 
positive for an order parameter with s-wave symmetry and zero for d-wave symmetry. We conclude that in 
the BCS approximation superconductivity with d-wave symmetry (or with other similar order parameters) 
is neither promoted nor suppressed by the on-site repulsion. Therefore we have to go beyond the mean-field 
level to find out whether a superconducting ground state is conceivable in the repulsive Hubbard model. 
This is a very difficult problem since it necessitates a reliable calculation of the correlation energy. 

The question of superconductivity in the repulsive Hubbard model (and the related t-J model) has been 
attacked by several different methods, the perturbative Renormalization Group (RG) J6] |7), numerical 
calculations of the renormalized scattering vertex (§1 and variational wave functions It9l ITOl ITTl fT2l 1131 . 
All these approaches are consistent with d-wave superconductivity in some doping range, but most of them 
have a very limited range of applicability. Thus the perturbative RG methods (6][3 are only valid for very 
small values of U and cannot provide a quantitative answer for the energy gap or the condensation energy. 
Conversely, studies of the t-J model ||9l[l0l[l2j[l2l can only be applied to the large U limit of the Hubbard 
model. Unfortunately, reliable results for the crossover region of the two-dimensional Hubbard model are 
still scarce. 



4 A refined variational wave function In our own work lfT4l we use a variational ansatz, which is a 
combination of Eqs. (f5]l and (O and constructed in such a way as to approach the crossover region from 
the Fermi liquid side and to allow for a superconducting ground state with d-wave symmetry, 



|$) = e-' lf e- 9£> \dBCS) . (9) 

Here the double occupancy is first partially suppressed (variational parameter g), and subsequently both 
the hopping of holes and the kinetic exchange are enhanced (variational parameter h). The additional 
parameter h reduces significantly the total energy, in agreement with similar earlier work |fT31 [161 . The 
parent state |g?BCS) is a BCS state with c?-wave symmetry, i.e., 

u k = -(H — ), u k v k - 



2 V ' E k h K K ~ 2£ k ' 



with 



Ek = J (ek - A 1 ) 2 + A k = A • (cos k x - cos k y ). 



Two further variational parameters are included in |dBCS), the BCS gap A and a "chemical potential" /i 
that allows to fix the average number of electrons. We emphasize that /i is not identical to the true chemical 
potential. 
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Fig. 1 (a): Total energy per site of the Gutzwiller-type wave function, (b): Total energy per site of the refined 
variational wave function. In both figures the lower curve corresponds to a precisely fixed number of electrons, while 
the upper curve represents the "grand canonical" case with fi tuned to represent a fixed average density. The calculation 
was carried out on a 8x8 square lattice for n = 0.8125 and U = St. 



The variational parameters are determined by minimizing the total energy (expectation value of the 
Hamiltonian). The natural way to get this energy for a given density is to project the parent state on a fixed 
number of electrons, working either in real space or in momentum space, and then to compute the energy 
using a Monte Carlo simulation. For the case h — (Gutzwiller type wave function), the parent state is 
easily written as a superposition of configurations in real space, where the number of doubly occupied sites 
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D is diagonal. The total energy, given by the lower curve in Fig. OJa), has a pronounced minimum around 
0.05<. The condensation energy is of the order of O.OOli. 

The case h > is more demanding because the operators T and D are diagonal with respect to dif- 
ferent bases. In a first approach, we have written the BCS state as a superposition of momentum space 
configurations with a fixed number of electrons. A discrete Hubbard-Stratonovich transformation is used 
to decouple the fermion operators in the Gutzwiller projector. The price to pay is a summation over "Ising 
spin" configurations, in addition to the summation over the momentum space configurations. This leads to 
a weight which is no longer positive definite in the Monte Carlo simulation. The resulting sign problem 
becomes severe as the gap becomes large. The result is represented by the lower curve of Fig. 0j6). In 
the range where our results are valid (do not suffer sign problems), the energy first decreases, reaches a 
minimum at about 0.03i and then increases. For A = O.OQt the average value of the sign is about 0.3. 
Results for larger values of A cannot be trusted. 

An alternative approach is to start with an unprojected parent state (i.e., with a fixed "chemical potential" 
instead of a fixed number of electrons). In this case only the summation over "Ising spin" configurations has 
to be performed, but the "chemical potential" has to be fixed to get the right average density. Interestingly, 
working with this wave function allows to escape the minus sign problem. The resulting total energy is 
represented by the upper curve of Fig. QJ&). The qualitative behaviour of the energy is quite different for 
the unprojected wave function; the minimum energy is located at A = 0! This result is confirmed by the 
same calculation carried out for h = (see Fig. QJa), upper curve). 

The two approaches should be equivalent in the thermodynamic limit. Unfortunately, for the system 
size studied so far (8 x 8), finite size effects seem to be rather large and, especially, lead to two conflicitng 
results. Larger system sizes together with finite size scaling will be required before we can reach any 
definite conclusion about the fundamental issue of superconductivity in the repulsive Hubbard model. 
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